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We solve Einstein's field equations coupled to relativistic hydrodynamics in full 3+1 general 
relativity to evolve astrophysical systems characterized by strong gravitational fields. We model 
rotating, collapsing and binary stars by idealized polytropic equations of state, with neutron stars 
as the main application. Our scheme is based on the BSSN formulation of the field equations. We 
assume adiabatic fiow, but allow for the formation of shocks. We determine the appearance of black 
holes by means of an apparent horizon finder. We introduce several new techniques for integrating 
the coupled Einstein-hydrodynamics system. For example, we choose our fluid variables so that 
they can be evolved without employing an artificial atmosphere. We also demonstrate the utility 
of working in a rotating coordinate system for some problems. We use rotating stars to experiment 
with several gauge choices for the lapse function and shift vector, and find some choices to be 
superior to others. We demonstrate the ability of our code to follow a rotating star that collapses 
from large radius to a black hole. Finally, we exploit rotating coordinates to evolve a corotating 
binary neutron star system in a quasi-equilibrium circular orbit for more than two orbital periods. 

PACS numbers: 04.30.Db, 04.25. Dm, 97.80.Fk 



I. INTRODUCTION 

With the availability of unprecedented observational 
data, the physics of compact object is entering a partic- 
ularly exciting phase. New instruments, including X-ray 
and 7-ray satellites and neutrino observatories, are de- 
tecting signals from highly relativistic events in regions of 
strong gravitational fields around neutron stars and black 
holes. A new generation of gravitational wave interferom- 
eters is promising to open a completely new window for 
the observation of compact objects. The ground-based 
gravity wave observatories LIGO and TAMA are already 
operational and are collecting data, GEO and VIRGO 
will be completed soon, and a space-based interferome- 
ter LISA is currently under design. 

Given the small signal-to-noise ratio in these new 
gravitational wave detectors, theoretical models of likely 
sources are needed for the positive identification of the 
signal as well as for its physical interpretation ||l|. One 
promising technique for the identification of signals in 
the noise output of the detector is matched filtering, 
which requires accurate theoretical gravitational wave 
templates 1^. The need for such templates has driven 
a surge of interest in developing reliable techniques ca- 
pable of their construction. 

Compact binaries, i.e. binaries consisting of either 
black holes or neutron stars, are among the most promis- 
ing sources of gravitational radiation. Much progress has 
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been made in refining post-Newtonian point-mass ap- 
proximations. These are suitable for large binary separa- 
tions for which relativistic effects are sufficiently small 
and any internal structure can be neglected 0. At 
small binary separations, the most promising technique 
for modeling the inspiral, coalescence and merger is nu- 
merical relativity. 

Several other observed phenomena involving compact 
objects require numerical relativity for their modeling. 
One such example is Gamma Ray Bursts (GRBs). While 
it is not yet known what the origin of GRBs is, the central 
source is almost certainly a compact object Most 
scenarios involve a rotating black hole surrounded by a 
massive magnetized disk, formed by a supernova, or the 
coalescence of binary neutron stars To confirm or 
refute any GRB scenario requires numerical studies in 
full 3+1 relativistic magnetohydrodynamics. 

Another astrophysical scenario requiring numerical 
treatment is the formation of supermassive black holes 
(SMBHs). Among the scenarios proposed to explain 
SMBH formation are the collapse of a relativistic clus- 
ter of coUisionless matter, like a relativistic star cluster 
1^ or self-interacting dark matter halo or the col- 
lapse of a supermassive star Depending on the de- 
tails of the collapse, SMBH formation may generate a 
strong gravitational wave signal in the frequency band 
of the proposed space-based laser interferometer LISA. 
Understanding the SMBH formation route may shed key 
insight into structure and galaxy formation in the early 
universe. 

Solving the coupled Einstein field and hydrodynamics 
equations is a challenging computational task, requiring 
the simultaneous solution of a large number of coupled 
nonlinear partial differential equations. In addition to 
all of the usual problems of numerical hydrodynamics - 
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handling advection, shock discontinuities, etc - one en- 
counters the problems inherent to numerical relativity. 
The latter include identifying a suitable formulation of 
Einstein's field equations, enforcing a well-behaved co- 
ordinate system, and, if black holes are formed, dealing 
with spacetime singularities. 

The construction of self-consistent numerical solutions 
to the coupled equations of relativistic hydrodynamics 
and gravitation dates back to the pioneering work of May 
and White in spherical symmetry j|] (see also jl^ for 
a review). In one of the first attempts to perform nu- 
merical integrations in three spatial dimensions, Wilson, 
Mathews, and Marronetti |jll|, ^ (see |l|| for 
later corrections) tackled the binary neutron star prob- 
lem. They simplified Einstein's field equations by as- 
suming that the spatial metric remains conformally flat 
at all times. Their implementation of relativistic hydro- 
dynamics was based on earlier work by Wilson [ p6| and 
used upwind differencing to handle advection and arti- 
ficial viscosity to capture shocks. The first fully self- 
consistent relativistic hydrodynamics code, which treats 
the gravitational fields without approximation, was de- 
veloped by Shibata This code, based on earher work 
by Shibata and Nakamura adopts a Van Leer hy- 
drodynamics scheme [ p^ |2(| and also employs artificial 
viscosity for shocks. This code has been used in vari- 
ous astrophysical applications, including the coalescence 
and merger of binary neutron stars |2^ and the sta- 
bility of single, rotating neutron stars p^ , p^ , [2^ . In 
an alternative approach, Font et al |2^ implemented a 
more accurate high-resolution shock-capturing technique 
to solve the equations of relativistic hydrodynamics. This 
code has been used to study pulsations of relativistic stars 



rately stable differentially rotating stars in equilibrium. 
We also show that our code can follow the collapse of 
rapidly differentially rotating stars reliably until an ap- 
parent horizon appears, by which time the equatorial ra- 
dius has decreased from its initial value by more than a 
factor of ten. 

We then turn to simulations of binary neutron stars. 
We adopt initial data describing corotating n — 1 poly- 
tropes in quasi-equilibrium circular orbit, and evolve 
these data for over two orbital periods. In this paper 
we present results for one particular binary and discuss 
the effect of corotating frames as well as the outer bound- 
aries. An extended study, including binary sequences up 
to the dynamically identified innermost stable circular 
orbit (ISCO), will be presented in a forthcoming paper 



This paper is organized as follows. Sees. || and III de- 
scribe our method of evolving the field and hydrodynamic 
equations, respectively. Sec. [V summarizes the various 
gauge choices with which we experiment. Sec. ^ lists the 
diagnostics used to gauge the reliability of our simula- 
tions. Sec . VII describes several tests of our algorithm. 



Sec. VIII applies our formalism to evolve non-rotating. 



uniformly rotating, and differentially rotating polytropes. 
In Sec. [X sketches our binary neutron star calculations. 
Our results are summarized in Sec. Some details of 
our hydrodynamic scheme and the rotating frame formal- 
ism are presented in the appendices. 



II. GRAVITATIONAL FIELD EVOLUTION 



A. Basic Equations 



In this paper, we report on the status and some astro- 
physical applications of our new 3-f 1 general relativistic 
hydrodynamics code. Our code, based on the so-called 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formula- 
tion of Einstein's equations has several novel 
features, including an algorithm that does not require 
the addition of a tenuous, pervasive atmosphere that is 
commonly used in Eulerian hydrodynamical codes, both 
Newtonian and relativistic. This "no atmosphere" al- 
gorithm proves to be very robust and eliminates many 
problems associated with the traditional atmospheric ap- 
proach pgt . 

We treat ID shocks, spherical dust collapse to black 
holes, and relativistic spherical equilibrium stars to 
demonstrate the ability of our code to accurately evolve 
the coupled field and hydrodynamic equations in rela- 
tivistic scenarios. We then use the evolution of stable 
and unstable uniformly rotating polytropes as a testbed 
to determine which gauge conditions are best-behaved in 
the presence of strong-field matter sources with signifi- 
cant angular momentum. We introduce rotating coordi- 
nate systems and show that these can yield more accu- 
rate simulations of rotating objects than inertial frames. 
We demonstrate the ability of our code to hold accu- 



We write the metric in the form 

ds^ ^ -a^de + {dx' + 13' dt) (dx^ + P^dt), (1) 

where a, and jij are the lapse, shift, and spatial met- 
ric, respectively. The extrinsic curvature Kij is defined 
by 



-2aKi. 



(2) 



where £p is the Lie derivative with respect to We 
choose geometrized units with G = c = 1 throughout, so 
Einstein's field equations are 



(3) 



We use greek letters to denote spacetime indices, and 
latin letters for spatial indices. Using the above vari- 
ables, the field equations (||) split into the usual 3-1-1 
ADM equations These consist of the Hamiltonian 

constraint 



R~ K,jK^^ + = 167rp, 



the momentum constraint 



D.Kl - D,K = %TTSi, 



(4) 



(5) 
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and the evolution equation for Ki. 
{dtKij - CpKij) = -DiDja 



+a{Ri 



2K,,K 



1 



KK,j (6) 
7.,(p-5%))) 



Here D, Rij and R are the covariant 



in addition to ( 
derivative operator, the three-dimensional Ricci tensor 
and the scalar curvature associated with jij . The matter 
source terms p, Si, and Sij are projections of the stress- 
energy tensor with respect to the unit normal n" on the 
time slice 



Sij = JiaJjpT"^. 



(7) 



Since numerical implementations of the ADM equations 
typically develop instabilities after very short times, we 
use a reformulation of these equations that is now often 
referred to as the BSSN formulation [|l8[ This re- 
formulation consists of evolving the conformally related 
metric 7^ , the conformal exponent (f), the trace of the 
extrinsic curvature K , the conformal traceless extrinsic 
curvature Aij , and the conformal connection functions 
defined by 



1^J 



(8) 
(9) 
(10) 



where det(7y ) = 1 and tr(^ij) = 0. In terms of these 
variables, Eqs. (0) and (ra) become 



{dt - Cfj)ji-i 
{dt - C0)<j> 



-2a A, 



= aK 
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idt-Cf3)K = -fW,D,a + -aK' 

+aA,jA'-^ + Ana{p + S) 
{dt-C0)A, = e-^^i-D,D,a + a{R,, -8t:S,j)) 



and 



-a{KAj-2AaA'j) 



(11) 
(12) 

(13) 

TF 

(14) 



(15) 



2~ 
'3 



2A'^dja 



(see 



6A'^0j - r^kA^'' + Snf^Sj^ 
for the computation of the Lie derivatives.) 



-2a ( - f'K^j 



In terms of the BSSN variables, the constraint equa- 
tions d) and (|) become, respectively. 



= H = fWiDje'f 



o50 



(16) 



^ 12 



where S"* = j"^^ Sj. While the two constraints are iden- 
tically zero for analytical solutions, they vanish only ap- 
proximately in numerical calculations. Thus, the Hamil- 
tonian and momentum constraint residuals Ti and M can 
be monitored as a code test during numerical evolution 
calculations. In the BSSN formulation, we also monitor 
the new constraint 



(18) 



B. Boundary Conditions 



Like any other hyperbolic system, the Einstein field 
equations must be supplemented by initial conditions 
and boundary conditions to have a unique evolution. 
We adopt boundary conditions that follow from the as- 
sumption of asymptotic flatness, i.e. gap rjap- In 
the asymptotic domain, monopole terms dominate in the 
longitudinal variables, so oc . The transverse fields 
will be dominated by outgoing gravitational waves, so 
lij — rjij (X f{t — r)r~^ and Aij (x g{t — r)r~^ , where / 
and g are unknown functions of retarded time. Note that 



is a special case of f{t 



so that (p, %j , and 



Aij all satisfy outgoing wave boundary conditions. The 
appropriate boundary conditions for K and depend 
on the gauge conditions used in the interior. 



C. Numerical Implementation 

We evolve Eqs. (pT|)-(p^) using an iterative Crank- 
Nicholson scheme with one predictor step and two cor- 
rector steps 1^ . In this algorithm a function / with time 
derivative / is updated from its value /" at timestep n 
to its value /"^^ at the next timestep rt + 1 a time AT 
later. In the explicit predictor step = /" + AT/", 

where /" is computed from quantities on timestep n, a 
"predicted" new value is found. In the following 

two corrector steps, V""^^ = /"-)-AT(/"-hi/"+i)/2 and 
jn+i ^ 3fn+i ^ /» + AT(/» + 2y«+i)/2, thcsc predicted 
values are "corrected". The final value /"^^ converges 
quadratically in AT. AT is set by the Courant factor: 
C — AT/ Ax, where Ax is the coordinate distance be- 
tween adjacent gridpoints. We typically use C = 0.5. 
The code implementing this evolution scheme has been 
discussed elsewhere |28|, so we will highlight here only 
the new features of our code. 
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We enforce the algebraic constraints det ^ij = 1 and 



tr(Ay ) = as described in 



Also following l33l, we 



replace the term ^f'^P^j in Eq. ( [Tsl ) with the analytically 
equivalent — (7*-' j + ^r*)/3^ j . These changes have little 
effect on the evolutions described in this paper, but lead 
to significant improvements when treating black holes by 
excision boundary conditions [^3| . 

We use second order centered differencing for all spa- 
tial derivatives in the field equations. We have not found 
it necessary to use upwind differencing for any deriva- 
tives. We did find, however, that the addition of some 
dissipation in the evolution equation for (j> increases the 
stability of the code. This can be supplied by upwind 
differencing of the term which advects cj) along the shift, 
and we have confirmed that this will indeed improve the 
stability. However, we have chosen instead to add the 
Hamiltonian constraint to the evolution equation for (p, 
as follows 



-—aK 
6 



chH. 



(19) 



Here the parameter ch is set between 0.02AT and 
0.06AT. cnTi. is a diffusive term, with Courant con- 
dition given Q by 2cif AT/(Aa;)^ < 1, so making ch 
proportional to AT is necessary in order to avoid an in- 
stability at high resolutions. It also provides dimensional 
consistency in (|l^). Using Eq. (|l^) offers the advantage 
of significantly decreasing the growth in the error of the 
Hamiltonian constraint (see Sec. [X for an example of 
this.) We note that the above is similar to one of the 
modifications of BSSN suggested in [p5|. 



D. Implementation of Boundary Conditions 



As discussed in Section ( |IIB| ), we use Sommerfeld 
boundary conditions for most of the field variables. That 
is, the value of a quantity / on the boundary at time t 
and distance r from the origin is 



/(r- Ar,t-AT), 



(20) 



where AT is the timestep and Ar = ae~^'^AT. 

For the functions F' we have experimented with sev- 
eral boundary conditions. We find little sensitivity to 
the condition used; the best choice seems to be fixing F* 
at their initial values (zero, for most of the applications 
here) . 



III. RELATIVISTIC HYDRODYNAMICS 

A. Basic Equations 

We describe the matter source of the Einstein equa- 
tions as a perfect fluid so that the stress-energy tensor 



can be written 



(21) 



Here poi £1 a-nd are the rest-mass density, speciflc 
internal energy, pressure, and fluid four-velocity, respec- 
tively. We adopt a F-law equation of state 



P= (F-l)poe, 



(22) 



where F is a constant. For isentropic flow, this is equiv- 
alent to the polytropic relation 



P 



(23) 



where k is a constant. In our simulations we encounter 
non-isentropic flow (due to shocks), and hence use equa- 
tion (|2|). 

The equations of motion follow from the continuity 
equation 

V^(pou^) = (24) 

and the conservation of stress-energy 

T^"",^ = . (25) 

Following ||l^, these equations can be brought into the 
form 



dtPi, + d,{pi,v') = 
dtCi, + i9^(e*w') = 

+ 2wh ^ 

(P,k , 

w 



(26) 
(27) 
(28) 



poau^e^''', w 



where h = 1 + e + P/ po, p* 
e* = {poe)^^^au'^e^''', Sk = pi,huk, and v 



-- p->,au", 
/vP is 

the 3-velocity. The quantity w is determined by the nor- 
malization condition w'^u^ = —1, which can be written 



w~ 



p^iwe^'l' / p^) 



r-i 



(29) 



The perfect fluid given by Eq. (|21|) generates the fol- 
lowing source terms for the ADM equations 



P 



w 



. SiSj 

h 



Pl^: 



(30) 
(31) 

(32) 



We will only be considering systems where there is vac- 
uum everywhere outside the star or stars. Therefore, the 
appropriate boundary condition on the matter flow is 
that no material should be flowing into the grid through 
the outer boundaries. 



5 



B. Numerical Implementation 



We evolve the hydrodynamic variables using an iter- 
ative Crank-Nicholson scheme. This scheme is slightly 
different from the one used to update the field variables. 
In the corrector steps, instead of weighting /" and 
equally (i.e. = p + AT(0.5/" + 0.5 7"+^)), we 

make the evolution more implicit by setting = 
/" + Ar(0.4/"-f 0.6 This makes the code slightly 

more stable. 

As is often done in hydrodynamics codes [^6| , the up- 
dating of the fluid variables onto a new timestep is di- 
vided into two steps ( "operator splitting" ) : the advection 
step (accounting for the advective terms on the left-hand 
sides of Eqs. (|2^)-(|2^)), and the source step (accounting 
for the right-hand sides of Eqs. (p6|)-(p^)). Each step of a 
Crank-Nicholson update consists of applying first an ad- 
vection substep and then a source substep. Our scheme 
for carrying out the advection substep is similar to the 
Van Leer scheme, and is discussed in detail in Appendix 
^. Since Eq. (|2^) has no sources, is completely up- 
dated after it is advected. Following we then use 
the updated to complete the updating of and Sk- 
it is shown in [ p9| that this gives improved behavior in 
Newtonian simulations of binary polytropes. 

C. Artificial Viscosity 
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FIG. 1: Maximum rest mass density po as a fraction of its 
initial value poi for the binary system shown in figure 
Stellar radial oscillations can be efficiently quenched by the 
proper use of linear viscosity, as shown here. The solid line 
shows the evolution without linear artificial viscosity, while 
the dashed line shows the efi'ect of this dissipative term. 

deep inside the neutrons stars. The (small) dissipated 
kinetic energy goes into thermal energy. We typically 
use 0.1 < Cqvis < 1-0. Linear artificial viscosity is not 
used in the runs described below. 



We handle shocks by adding quadratic artificial viscos- 
ity. This consists of adding a viscous pressure 

p _ / CQ^isAiSvf for Sv<0 , , . 

Q^''' ~ \ otherwise , ^'^'^^ 

r 

where A is defined as ^^^6^/*^)r-i and 6v = 2dkv''Ax. 
Shock heating causes an increase in the local internal 
energy. Following Q, we change equation (|2^) to 

We have also implemented linear artificial viscosity 
terms that can be used to dissipate radial oscillations 
triggered in stars by the truncation error associated with 
finite differencing. The corresponding addition to the 
pressure is 

p ^ f -C-Lyis\/{T/n)p^A 6v for 5v < , 
\ otherwise . 

Linear viscosity can be used at the beginning of a run 
to drive the initial data to dynamical equilibrium and 
later switched off. Figure ^ shows an example of how 
the radial oscillations can be quenched by linear viscos- 
ity. For this particular example, the PlvIs was active 
only where the rest mass density exceeded a particu- 
lar threshold value, to force this dissipative effect only 



D. Non- Atmospheric Hydrodynamics 

Numerical work in Eulerian hydrodynamics, both 
Newtonian and relativistic, has typically required the 
presence of a pervasive tenuous "atmosphere" that covers 
the computational grid outside the stars. To our knowl- 
edge, most published codes to date need to keep a min- 
imum nonzero density that is usually set to be several 
orders of magnitude smaller than the maximum stellar 
density. Such an atmosphere has been necessary to pre- 
vent overflows arising from dividing by density in cells 
devoid of matter. This artificial atmosphere has to be 
small enough not to affect the true dynamical behavior 
of the system. However, very small values will propagate 
the round-off numerical error very quickly every time a 
division by the density is performed. A problem with 
the presence of this atmosphere is that as soon as the 
time evolution starts the material begins to fall onto the 
star, creating accretion shocks. Swesty et al. solved 
this problem by adding a non-zero temperature to the at- 
mosphere to restore some sort of equilibrium that would 
counterbalance the infall. Also, in order to avoid the 
bow shocks generated in the atmosphere by two stars in 
circular orbital motion, these authors provide the atmo- 
sphere with initial angular velocity. These are some of 
the typical problems present in the traditional artificial 
atmosphere approach found in many Eulerian hydrody- 
namics schemes. 
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In this paper we present a very simple algorithm that 
does not require the presence of atmospheric material. It 
consists of two ingredients. The first is the use of the 
spatial components of the linear momentum variable Sk 
as our hydrodynamical variable instead of the tra- 
ditional fluid four-velocity spatial components used in 
most hydrodynamical codes (see for instance (l^). In 
the latter case, the Euler equation is used to update the 
flux {pi, M*). Once this update is completed, the dynam- 
ical field It* is recovered by dividing by the density p^. 
Using Sk as a variable, we avoid these divisions. The 
only time when the variable needs to be calculated 
explicitly is when we need the three-velocity that ap- 
pears in every advection term on the left hand side of 
Eqs. (HIH). To avoid doing this calculation for very 
low values p*, we add the second ingredient: the intro- 



at the outer-x boundary i 



is implemented as 



duction of a threshold value 



below which all the 



hydrodynamical fields are set to vacuum values (i. e. 
Pi, = = = 0). A typical value for p^, is 10"'' 
times the maximum initial value of Pi,. 

However, as the time evolution progresses, a tenuous 
shell of material typically drifts away from the stars and 
creates regions of very low density outside our stars. If 
nothing special is done about them, small shocks will 
heat this low-density region to very high temperatures 
and will generate large velocities. Although this low- 
density region has a negligible effect on our stars and 
spacetimes, it can cause the code to crash. Therefore, we 
impose a heating limit outside the star 

= min(e^, lOp*) if < efactor x p^max, (36) 

where efactor is a constant that is determined empirically 
for a given physical scenario. We generally choose values 
between 10~^ and 10~^, where the larger values of efactor 
were only needed in simulations of collapsing stars with 
a strong bounce. We note that this is similar to the 
technique used in [p6| , in which the polytropic equation 
of state is applied in the low-density region outside 
the star or stars. 



E. Boundary Conditions 

Since matter often diffuses outward, albeit in minute 
quantities, from the surface of the star(s) to the bound- 
aries, we need to impose boundary conditions on the 
matter at the outer grid points. In algorithms where 
an artificial atmosphere is present, it is crucial to choose 
boundary conditions which do not lead to a continuous 
inflow from the boundary, or to bad behavior in the at- 
mosphere. By eliminating such an atmosphere, however, 
all reasonable boundary conditions yield the same behav- 
ior so long as the boundaries are placed far enough from 
the star(s) that little matter ever reaches them. 

We usually use an outflow boundary condition. For ex- 
ample, if the x-coordinate of gridpoints is indexed by an 
integer i with imi-a < * ^ *max, this boundary condition 



n+l 
★imax 
n+l 
★imax 



— / "^irnax-lif '^'j^^^-l ^ 







otherwise 



(37) 
(38) 

(39) 



We have experimented with other boundary conditions 
as well. We have tried fixing p^,, e^,, and Sk at their initial 
values. We have also tried simply copying the adjacent 
gridpoint onto the boundary with no outflow restrictions 
{Copy). These conditions produce similar results to those 
of the outflow condition for all applications, while being 
somewhat less computationally expensive. 



IV. GAUGE CHOICES 
A. Lapse 

We experiment with several time slicing conditions. 
First, we try maximal slicing, which enforces K = dtK — 




: 



-j'WjD.a + aA.jA'^ + Awaip + S) . (40) 



This slicing condition has the advantages of controlling K 
and avoiding singularities. Unfortunately, it is a compu- 
tationally expensive gauge choice, since it involves solv- 
ing an elliptic PDE every timestep. Therefore, we also 
try a slicing condition which approximates maximal slic- 
ing, the so-called "K-driver" proposed by Balakrishna et 
al ]39|] . The idea is to convert the elliptic equation (max- 
imal slicing) into a parabolic evolution equation 



Kdr: dta = -e{dtK + cK) , 



(41) 



where e and c are positive constants. The equation 
dtK = —cK, corresponding to exponential decay in K, is 
the solution of equation (|4^) as e — > oo. However, setting 
e at too large a value in our code will produce a numerical 
instability. (See the discussion of ch in Sec. 11 C.) For- 
tunately, this limitation can be overcome. We are able 
to effectively evolve with larger e by breaking up each 
timestep into several substeps and evolve Eq.(^l|) using 
a smaller AT than that used by the other variables. On 
each substep, we use the values of the metric on the des- 
tination time level, so the process is equivalent to solving 
the elliptic equation dtK + cK = by relaxation, except 
that we do not carry the process to convergence. Instead, 
we typically use 5 substeps per step, with e = 0.625 and 
c = 0.1. 

An even less computationally expensive lapse condition 
is harmonic slicing, which for vanishing shift reduces to 



hm: (9, 



-1^-1/2) 







(42) 



We apply this condition unchanged for vanishing and 
non-zero shift, and find that it often gives behavior sim- 
ilar to that obtained by using the above two slicings. 
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B. Shift 

We also experiment with different spatial gauge 
choices. The simplest admissible shift choice, which turns 
out to be surprisingly good for collapsing star applica- 
tions, is to keep the shift "frozen" at its initial values 



fz: p\t)^p\0) 



(43) 



at each grid point. 

We also try the approximate minimal distortion 
(AMD) gauge introduced by Shibata |40| 



AMD: %V2/3» + = J, 



(44) 



where is the flat-space Laplacian and 



J,; = 167raS', + 2A,j{a^^ - 6a0'^) -I- -aK,^. (45) 

o 

This gauge condition was designed to approximate the 
Smarr and York minimal distortion shift condition 
which in turn was constructed to minimize gauge-related 
time variation in the spatial metric. 

As Shibata points out |Q, the AMD condition must 
be modified in the event of a collapse in order to prevent 
a "blowing out" of coordinates on the black hole throat, 
which manifests itself by a growth in the proper 3- volume 
element, i.e. by growth in 0. The blowing out can be 
controlled by preventing the radial component of the shift 
from becoming large and positive 



MAMD: = 




for < (4/3)(/., 



//3amd— otherwise 
r 



(46) 

where is the solution of equation (|4J), (p^ is the 

value of (j) at the coordinate origin, 0ci = (f'cit = 0), and 



/ 



30c 



1 



1 + {r/RY 



(47) 
(48) 



where i? is a constant. This correction is only useful in 
configurations with near spherical symmetry, so that the 
collapse is nearly radial at the center. It is disabled for 
simulations of binary systems. 

Finally, we try approximating the "Gamma-freezing" 
condition dtT^ — using a "Gamma-driver" , which en- 
forces controls in the same way that the K-driver con- 
trols K 



Gdr: 9*/?* = /fc(9tr + ryP ). 



(49) 



Here k and rj are positive constants, and, as with the 
K-driver, we can effectively make k larger than would 
otherwise be possible by breaking up each step into mul- 
tiple substeps. This shift condition has been used suc- 
cessfully in black hole evolution calculations [E2|. The 



Gamma-freezing condition is closely related to minimal 
distortion (and hence approximate minimal distortion) , 
and it is hence not surprising that the modification ( |4q ) 
must also be applied to the Gamma-driver shift. Typical 
values for the Gamma-driver's parameters are k — 0.01 
and rj — 0.2, using 10 substeps per step. 

We have also implemented of the full Gamma freez- 
ing condition (i9tF' = 0). However, applying this con- 
dition requires solving three coupled elliptic equations 
each timestep (see Eq. (|l5|)), and we have found Gamma 
freezing to be too computationally expensive to be worth 
solving exactly. 



C. Rotating Frames 

Rotating coordinate frames possess superior angular 
momentum conservation capability over inertial frames 
in many applications, such as the hydrodynamical evolu- 
tion of binaries systems. In transforming from an inertial 
frame with coordinates (i, x, y, z) to a rotating frame with 
coordinates {t,x,y,z) and constant angular frequency 
il — flel, we apply the following relations 

t = t 

X — X cos(51t) — y sin(r2t) 

y = X sin(r2i) -I- y cos(f2t) 

z = z, (50) 

where the barred variables will represent quantities in 
the inertial frame in the remainder of this section. It is 
convenient to compare variables in the two frames at an 
instant t = t = Q a,i which the two frames are aligned. 
At this instant, the line element transforms from 

ds^ ~{a - (3'P,)dP + 2(3idx'dt + %dx'dx' 

to 

ds^ = ~(a-%0' + inxr)'){0i + {nxf)^)^dt^ 

+ 2%{P'' + (d X r)^)dxHt + %dx'dx3, 

where r= {x,y,z). From this equation, we sec that the 
following transformation rules apply at f = < = 0: 

a — a 

= p' + {nx fy 

7u = (51) 

Eq. (^ ) provides the transformation rules for the initial 
metric data from the inertial frame (where it is usually 
derived) to the rotating frame. The only change is the 
addition of a new term in the shift. At later times, vectors 
and tensors in the two frames will also differ by a rotation. 
However, we note that at all times there will be some 
inertial frame, related to (t, a;, y, z) by a rotation matrix, 
which has axes aligned with the rotating frame and whose 
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metric is related to that of the rotating frame by Eq. (51). 
Using the coordinate transformations (|5^) we can derive 
the relation between all the fields in both frames, for 
example 



u 

Ui 



u 



(f1 X r)\ 



(52) 



where is the fluid three- velocity = /vP . Ait — i — 
0, the components of any spatial tensor are unchanged 
under this transformation, for example 



7-y =7y, 



(53) 



and equivalently for and ify. The relation ( |5^ ) im- 
plies = 7*-' , since the inverse of a tensor is unique, so 
that we also find K^^ — K'^^ . To complete our introduc- 
tion of rotating frames in general relativity, we refer the 
reader to Appendix where we show that the Newto- 
nian limit of the relativistic Euler equation with a shift 
vector of the form of equation ( |5l| ) reduces to the fa- 
miliar Newtonian form of the Euler equation in rotating 
frames. In Appendix we show that the integrands 
used to evaluate M, Mq, and J in Eqs. (|6^), (64), and 
), respectively, remain unchanged when expressed in 



(I 

terms of rotating frame variables. 



D. Boundary Conditions 

We always choose initial data which satisfy maximal 
slicing (^) and gauge choices which approximately main- 
tain this slicing. Far from the source, equation ( ^ ) be- 
comes the Laplace equation, and its solution can be writ- 
ten as a sum of multipole moment fields. In the pres- 
ence of matter, the source will always have a nonzero 
monopole moment, so the asymptotic form of the lapse 
is 



1 cx r 



(54) 



All of our spatial gauge choices resemble one another, 
so we will just derive the shift boundary condition for the 
AMD shift (Ej ) , which is the easiest. The three compo- 
nents of Eq. (44) can be decoupled by decomposing /3* 
as in [Eol 



Ip. 



where x'^ are the Cartesian coordinates. Eq. (^ 
becomes 



JiX . 



(55) 
then 

(56) 
(57) 



To lowest order, Ji = pv^. We will be studying systems 
with azimuthal velocity fields, for which = 0, and 



hence rj = = 0. The lowest nonvanishing moment of 
Ji, from the monopole piece of p, is ^ = 1, m = ±1. We 
can solve the Laplace equation (outside the star) assum- 
ing asymptotic flatness to get the boundary conditions 
(3^ cx yr~^, fi^ oc xr~^, and, to this order, (3^ = 0. A 
nonzero boundary condition for (3^ must come from a 
higher-order term, but, since (3^ will be very small, our 
simulations are insensitive to it. We use 



f3^ oc yr 



„-3 



(3^ oc xyzr 



(58) 



The (3^ condition is obtained by ignoring 77 and solving 
( p6| ) subject to the lowest-order moment of the Aija'^ 
term in Ji (see Eq. (|5|).) 

Note that the coordinate-rotation component of a shift, 
i1 X r, is a homogeneous solution of equation ( ^ ) . It was 
eliminated in (^) by the assumption of asymptotic flat- 
ness. When working in a rotating frame, the shift does 
not obey an asymptotic flatness condition. Indeed, one 
can think of a coordinate rotation as a boundary condi- 
tion imposed on the shift. In such frames, the asymptotic 
form of the shift is l1 x r plus a piece which behaves like 
(p8|), and the boundary conditions must be set accord- 
ingly 



V. DIAGNOSTICS 

In order to gauge the accuracy of our simulations, we 
monitor the L2 norms of the violation in the constraint 
equations. These are the Hamiltonian constraint 7i (0), 
the momentum constraint Ai^ (pT|), and the Gamma con- 
straint (^) . We normalize the Hamiltonian and momen- 
tum constraint violation by their L2 norm by 



N 



HC 



(2^i^V)'+(^^A^)% (59) 



2\ 



and 



MC 



E 



Id'K ] (60) 



1/2 



The two terms in the Gamma constraint (|l^) often van- 
ish individually, so that a similar normalization is not 
meaningful for this constraint. 

Related to the Hamiltonian and momentum con- 
straints are mass and angular momentum conservation. 
In Gartesian coordinates, the ADM mass M and the an- 
gular momentum J* are defined by the behavior of the 
metric on a closed surface at asymptotically flat spatial 
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infinity 



M ^ -Ir f V77™7'"(7mnj-7jn,™)d'5.(61) 

(62) 



J,. = 



S„ 



Since Eq. (^2|) is computed at spatial infinity, e^j'^ is 
the flat-space Levi-Civita tensor. Using Gauss' Law, we 
transform the surface integrals into volume integrals 



iy V 167r ' 247r 



(63) 



„ 1 p^jTcp 1^^\^3^ 
167r 167r / 

"^^ = ^^i' J^ii^^i+^'Sk (64) 

(see, for example. Appendix A in [ ^Sf for a derivation). 
Note that since Ey'^ is outside the integral, it is still 
the flat-space Levi-Civita tensor. Baryon conservation 
(pom^);ai — implies that the rest mass 



Mo = / p^d^x 



(65) 



is also conserved. Due to the finite differencing in our 
hydrodynamic scheme, Mq is conserved identically except 
for matter flow off the computational grid. We therefore 
monitor AIq only as a diagnostic of how much matter 
flows through the outer boundaries. 

Eqs. (|6^) and ( |6^ ) are only valid in asymptotically flat 
spatial hypersurfaces and thus are not suited for use in 
rotating reference frames. However, the problem can be 
sidestepped quite easily by calculating the mass and an- 
gular momentum in the inertial frame, as functions of the 
dynamical variables of the rotating frame. In Appendix 
y we show that the integrals for these conserved quanti- 
ties are exactly the same when expressed in terms of the 
rotating frame fields. 

Another useful quantity to monitor is the circulation. 
According to the Kelvin-Helmholtz theorem |43j , the rel- 
ativistic circulation 



C(c) 



hUf^X^da 



(66) 



is conserved in isentropic flow along an arbitrary closed 
curve c when evaluated on hypersurfaces of constant 
proper time. Here h— 1 + e + P/pis the specific en- 
thalpy, fj is a parameter which labels points on c, and 
is the tangent vector to the curve c. Since C(c) is 
only conserved for isentropic flow, checking conservation 
of circulation along a few curves will measure the impor- 
tance of numerical and artificial viscosity on an evolution. 
We do not monitor circulation in this paper, although 
such a check has been implemented elsewhere |Q. 

Finally, we check for the existence of apparent horizons 
using the apparent horizon finder described in ES]. 



Proc 1 



Proc 3 



Proc 2 



Proc 4 



FIG. 2: This diagram shows how our code implements tt- 
symmetry in distributed-memory computer clusters. The 
black circles correspond to grid points, and the bottom row 
corresponds to the boundary in the plane orthogonal to the 
rotation axis. The white circles represent the ghostzones 
needed by our second order finite difference stencil. The ar- 
row connects two points that are related in the presence of 
TT-symmetry. 



VI. NUMERICAL CODE DESCRIPTION 

All our algorithms have been implemented in a paral- 
lel, distributed-memory environment using DAGH soft- 
ware developed as part of the Binary Black Hole 
Grand Challenge Alliance. When we need to solve ellip- 
tic equations (to construct initial data or to impose ellip- 
tic gauge conditions), we use the computational toolkit 
PETSc 0. 

Due to our large number of variables, the memory 
needed by our code is considerable (for example, a run 
with 64^ spatial zones may require up to 2 Gbytes of 
memory). Thus, it is crucial that we exploit any sym- 
metries present in a given problem to minimize the num- 
ber of gridpoints needed. We have implemented reflec- 
tion symmetry across a coordinate plane (equatorial sym- 
metry) and reflection symmetry across three coordinate 
planes (octant symmetry) , which cut the size of our grids 
by factors of two and eight. Our code also allows us to 
enforce 7r-symmetry, which assumes symmetry under a 
rotation of tt radians about a given axis. Unlike equa- 
torial and octant symmetry, the implementation of tt- 
symmetry is not trivial on distributed-memory parallel 
systems. This is because grid points needed to gener- 
ate the proper boundary conditions at a given location 
of the outer grid boundary will usually be located in the 
memory of a different processor, as seen in the diagram 
of Fig. 1^ where the value of the field at the white circle 
needed by point P in processor 4 must be provided by a 
black circle on processor number 3. We fix this problem 
by creating a two dimensional array for each field that 
stores the values of the field on all the grid points out- 
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side the boundary (white circles) needed to calculate the 
derivatives of the field at the grid points at the bound- 
aries (first row of black points). Each processor is re- 
sponsible for updating the array values corresponding to 
grid points within its domain by a 7r-radian rotation. Up- 
dated values are broadcast via MPI, and each processor 
has a copy of the complete two-dimensional array from 
which to draw the corresponding boundary values. 

VII. TESTS 

A. Vacuum Code Tests 

The algorithm for evolving the field equations was first 
tested in the context of small amplitude gravitational 
waves With harmonic slicing, the system could 

be accurately evolved for over 100 light crossing times 
without any sign of instability. The results also showed 
second-order convergence to the analytic solution when 
resolution was increased. In a forthcoming paper 
it will be demonstrated that this same code can stably 
evolve isolated black holes, with and without rotation, in 
Kerr-Schild coordinates. 




-1 -0.5 0.5 

X 



FIG. 3: The one-dimensional relativistic Riemann shock tube 
test. We plot the numerical rest density po (triangles), pres- 
sure P (squares), and velocity v (crosses) at t = 0.5. Solid 
curves show the analytic values. This particular run used 

Cqvis ~ 1. 



B. Hydro-without-Hydro 

Next, it was demonstrated that this field evolution 
scheme is stable when predetermined matter sources are 
present |48| . This was done by inserting the matter 
sources from known solutions of the Einstein equations 
and then evolving the gravitational field equations. Using 
this "hydro-without-hydro" approach, evolved the 
Oppenhcimer-Volkoff solution for static stars without en- 
countering any instability, and the Oppenheimer-Snyder 
solution for collapse of homogeneous dust spheres well 
past horizon formation. The same hydro-without-hydro 
approach was later used to model the quasi-equilibrium 
inspiral of binary neutron star systems and calculate the 
complete late-inspiral gravitational wavetrain outside the 

isco Bgl, |50|1. 



C. Shock Tube 

Every hydrodynamic algorithm must demonstrate 
some ability to handle shocks. In Fig. ||, we compare the 
output of our code for a simple one-dimensional shock 
tube problem with the exact result, which is known an- 
alytically in special relativity |pl| . In order to compare 
with this result, the metric functions are held at their 
Minkowski values throughout this test. At t = 0, we set 
V = — everywhere. For x < we set po = 15, 
P — 225 initially, and for a; > we set po = 1, P — 1. 
We output data at t = 0.5. In Figure 1, we use ar tificia l 
viscosity parameters Cqvis — 1, ClvIs — (see Sec. [II C] ) 
and a grid with 400 points. The shock is resolved quite 



well, and the only disturbing feature of our results is the 
"overshoot" in variables at the rarefaction wave. Nor- 
man and Winkler have shown that these overshoots 
are present in the solution to the finite difference equa- 
tions of artificial viscosity schemes even in the limit of 
the grid spacing going to zero. This problem therefore 
represents a fundamental limitation of artificial viscosity 
schemes, and points to the need for more sophisticated 
high- resolution-shock-capturing techniques when strong 
shocks are present (see, e.g., p6|). However, for many 
of our astrophysical applications (e.g. binary inspiral) we 
anticipate at most very weak shocks, so that the use of 
artificial viscosity schemes is adequate. 

Our results are completely insensitive to CqvIs when it 
is within the range — 0.1. We find the optimal behavior 
around Cqvis ~ 1, at which point the effects of artifi- 
cial viscosity are small but noticeable. For Cqvis ~ 5 or 
greater, the viscosity is too large, and we are unable to 
evolve accurately. 

Note that in the above example > lO^^p^max ev- 



erywhere, so e*-limiting ( p6D is never used. More ex- 
treme shocks can be created by increasing the density 
ratio po{x > 0)/po{x < 0). We find that we can treat 
shocks reasonably accurately for ratios of up to about 
twenty. 



D. Oppenheimer-Snyder dust collapse 

As a second simulation which can be tested against ex- 
act results, we model Oppenheimer-Snyder (OS) collapse 
of a homogeneous dust sphere to a black hole ||5^. The 
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FIG. 4: The evolution of the conformal exponent at the 
origin (pc and the ADM mass Af during Oppenheimer- Snyder 
collapse. The deviation of from its analytic value <^c"^' is 
measured by A(f)c = (0c - (fiT'^^) / + This is plotted 

on the top panel. On the bottom panel, we plot the ADM 
mass of the system and the irreducible mass of black hole, 
given by the area of the apparent horizon. We compare runs 
at two different resolutions. 



analytic solution for OS collapse can be transformed into 
maximal slicing and isotropic coordinates following [Q. 
We use the analytic solution at i = 0, when the mat- 
ter is at rest, as initial data for all variables. We then 
evolve the gravitational and hydrodynamic fields with 
our 3+1 code and compare the result with the exact so- 
lution. At each timestep, we determine the lapse by solv- 
ing the maximal slicing condition from the fields on our 
3D grid. For the shift we insert the analytic values cor- 
responding to isotropic coordinates. We evolve on a 32^^ 
grid and a 64"^ grid, utilizing octant symmetry to treat 
only the upper octant. Our outer boundaries are placed 
at 4Af in the isotropic coordinates of our grid. The initial 
Schwarzschild radius of the dust sphere is 3M. 

In Fig. ^, we show the convergence of the central con- 
formal exponent (j>c to the exact value. In Fig. ^ we com- 
pare the density profiles at several times for the 64^ grid 
to their analytical values. Throughout the evolution, we 
search for apparent horizons. At t = 8.75A/, we locate an 
apparent horizon with irreducible mass Mah/M — 1.03. 
(See Fig. ^.) This mass remains constant to within 3% 
until the end of the simulation (« 3Af later). As is known 
analytically, all of the mass falls inside the black hole. 

This test is similar to the "hydro-without-hydro" 
Oppenheimer-Snyder test performed on our code in [Q, 
except that here the matter fields and the lapse are de- 
termined numerically rather than set to their analytic 
values. 



FIG. 5: The density, defined by Equation (^) as a func- 
tion of isotropic radius during Oppenheimer-Snyder collapse. 
We compare our numerical results (crosses) with the analytic 
profiles. 



VIII. SINGLE STARS 

In this Section we study isolated stars, both non- 
rotating and rotating. The initial data, constructed from 
the OV solution for non-rotating equilibrium stars and 
with the code of [55| ] for equilibrium rotating stars, are 
summarized in Table || (see also Fig. |6l) We use the 



same coordinates as used in [55 (except transformed 
from spherical to Cartesian). For spherical, nonrotating 
systems (OV stars), these are the familiar isotropic coor- 
dinates. In these coordinates, the 3-metric for stationary 
spherical stars is conformally flat, and the event horizon 
of a Schwarzschild black hole is located at r = 0.5Af . All 
stars are n = 1, F = 2 polytropes (see Eq. (p3|)), and 
are dynamically evolved with the gamma-law equation 
of state (^). The nondimensional units throughout are 
set by requiring k — G = c — 1. 



A. Static Stars 

The stability properties of non-rotating F = 2 poly- 
tropes are known analytically and can be used as a test 
of our code. We use the OV |Q solution describing 
equilibrium polytropes in spherical symmetry as initial 
data, and evolve the matter and fields dynamically. An 
OV star is characterized by one parameter, which can be 
taken to be the central rest density pc- (We will hence- 
forth drop the subscript "0" on the rest density when 
referring to central rest density.) Along the sequence of 
increasing p^i the mass M takes a maximum value A^max 
at a critical central density p™*- (^^^ Fig- !§[) Stars 
with Pc < p"'* are dynamically stable, while stars with 
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TABLE I: Isolated Equilibrium Star Configurations (F = 2). 



Star 


PI" 


Mo" 




r> d 

Uc 




T/\WY 








A 


0.157 


0.171 


0.700 


0.866 


0.00 


0.000 




1.00 




B 


0.162 


0.178 


0.540 


0.714 


0.00 


0.000 




1.00 




C 


0.170 


0.186 


0.697 


0.881 


0.35 


0.032 


1.00 


0.88 




D 


0.171 


0.187 


0.596 


0.780 


0.34 


0.031 


1.00 


0.87 




E 


0.279 


0.304 


1.251 


1.613 


1.02 


0.230 


2.44 


0.30 




F 


0.049 


0.050 


1.240 


1.290 


0.72 


0.053 


5.88 


0.75 





ADM mass 

rest (baryonic) mass 

coordinate equatorial radius 

areal radius at the equator 

ratio of angular momentum to 

ratio of kinetic to gravitational potential energy 

ratio of central to equatorial angular velocity 

ratio of polar to equatorial coordinate radius 



0.1 



0.05 



-0.05 





1 1 1 1 1 1 1 1 1 1 1 1 




(16)3 




(32)3 




\ (64)3 







2 3 



0.18 



0.16 



0.14 



0.12 



(D) 



j=o.or-- 




0.2 



0.4 

Pc 



0.6 



FIG. 6: Stars A, B, C, and D and the constant-J sequences 
on which they lie. Open circles represent stable configura- 
tions, and closed circles denote unstable configurations. 



FIG. 7: Fractional change in the central rest density of star 
A when evolved on grids of three different resolutions. 




Pc > Pc^^* are unstable and collapse to black holes on a 
dynamical timescale. The dynamical timescale is given 

— 1/2 

by the free-fall time pc ■ To verify that our code can 
distinguish stable and unstable configurations we evolve 
two very similar models on either side of the critical point 
at p™*. 

In our units, p™' = 0.32 and A/,nax = 0.164. Star A 
has an initial central rest density pd = 0.2 and is there- 
fore stable. We set our outer boundaries at x,y,z = 2 
and evolve this star with three different resolutions 16"^, 
32'^, and 64'^, once again utilizing octant symmetry. In 
Fig. ^ we show the central density evolution for the three 
resolutions using harmonic slicing and the Gamma-driver 
shift. We see that our code does converge to the exact 
(stationary) solution. There are three sources of the de- 
viations from exact second-order convergence (see also 



FIG. 8: Star A evolved on a 32^ grid using various gauge 
choices. Here "hm" refers to harmonic lapse, "Kdr" to K- 
driver lapse, "fz" to frozen shift, and "Gdr" to Gamma-driver 
shift. 



p8[). First, there are components of the error which 
scale with a higher power of the grid width (e.g. Ax"^). 
Second, there is the noise caused by discontinuities at 
the surface of the star. Finally, errors are generated by 
imposing outer boundary conditions at finite distance. 

In Fig. ^, we evolve on a 32"^ grid with several gauge 
choices. Already we see that the choice of gauge is im- 
portant. Even in this static case, where the shift in the 
OV solution vanishes, it is necessary to use a dynamic 
shift for long term stability. With the Gamma-driver, we 
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evolve to tpc = 50 {t/M = 712) and never encounter 
an instability. 

We have stably evolved star A on the 32'' grid for 
many fundamental radial oscillation periods, which have 
a period of about = Ipc'^^'^ [0. However, we find 
that high-frequency, high-amplitude oscillations appear 
in after a few periods and persist thereafter. The 
onset of these oscillations can be delayed and their am- 
plitude diminished by increasing grid resolution. They 
can be removed altogether by making the hydrodynamic 
algorithm more implicit, i.e. by increasing the weight 
on th e new timestep in the corrector step (See Section 
[II B ). This adversely affects our ability to handle shocks, 
though. The problem may also be resolved by the use of 
a more sophisticated hydrodynamics scheme (see, e.g., 
). For the less relativistic stellar model used by 
27 1 our code produces non-physical high-frequency 
oscillations after about 6 radial oscillation periods r^. For 
the applications discussed here, these late-time problems 
are not relevant. 

Star B has an initial central density of pd — 0.4 and 
is dynamically unstable. We evolve this star with har- 
monic lapse and frozen shift, imposing outer boundaries 
at 1.2 M, on two different grids (32^ and 64^). In Fig|, 
we plot the central density and lapse as a function of 
time. The collapse is induced solely by the perturba- 
tions caused by putting the star on a discrete grid. Since 
these perturbations become smaller as grid resolution is 
increased, it is not surprising that the star on the lower- 
resolution grid collapses before the one on the higher- 
resolution grid. Since both collapse, it appears that 32'^ 
zones are sufficient to distinguish stable from unstable 
stars. Eventually, the star collapses to a point at which 
there are too few grid points across the star's diameter 
for the evolution to remain accurate. We terminate our 
evolutions when the error in the ADM mass exceeds 15% 
of the original mass. The 32'^ grid turns out to be too 
coarse for an apparent horizon to be located. We do lo- 
cate an apparent horizon in the 64"^ run shortly before 
the simulation is terminated. At this point the central 
lapse has collapsed to ac — 0.05, and, as a measure of 
error, the ADM mass deviates by 10% from its initial 
value. The horizon mass agrees well with the ADM mass 
M ^Mah- 

Also included in Fig. ^ is a 64^ simulation using har- 
monic lapse and the Gamma-driver. Similar behavior is 
seen in these coordinates. We will investigate the per- 
formance of various gauge choices in more detail in the 
following Section. 



B. Uniformly Rotating Stars 

1. Inertial Frame 

Simulations of systems with non-zero angular momen- 
tum are very sensitive to the choice of coordinates, which 
makes them very good test cases for comparing the nu- 
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FIG. 9: The collapse of star B seen with various gauge 
choices. The abbreviations are the same as in Fig. 
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FIG. 10: Star C evolved on a 64 x 32^ grid with Gamma- 
driver shift condition and various lapse choices. 



merical behavior of different gauges and slicings. Most 
of these effects can be seen when we evolve uniformly 
rotating stars. 

We consider two uniformly rotating stars, stars C 
and D, on one constant angular momentum sequence 
J — 0.01 (see Fig. The J = 0.01 sequence has a 
turning-point at - 0-4, M^ax = 0.172. For a se- 
quence of uniformly rotating stars, this turning point 
marks the onset of secular, not dynamical, radial instabil- 
ity |57| . It is possible for a star on the secularly unstable 
branch to be stabilized temporarily if the star begins to 
rotate differentially, so that no instability will develop on 
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FIG. 11: Star C evolved on a 64 x 32^ grid with K-driver 
lapse and various shift choices. We plot the maximum value 
of the absolute value of on the grid to show the dependence 
of F' on spatial gauge. 



a dynamical timescale. However, prior numerical stud- 
ies have found the point of onset of dynamical in- 
stability to be very close to the point of onset of secular 
instability, which we confirm with our simulations here. 

We again pick two similar stars on either side of the 
onset of secular instability: star C with pd = 0.3 on the 
stable branch and Star D with pd = 0.5 on the unsta- 
ble branch. We dynamically evolve these two stars with 
different choices for the slicing and gauge. All simula- 
tions are performed on 64 x 32^ grids, utilizing equato- 
rial and TT-symmetry to evolve only half of a hemisphere. 
The outer boundaries are placed at [—1.5, 1.5] x [0, 1.5]^. 
There are now two relevant timescales - the free-fall time 

— 1/2 

r// ~ /9c and the orbital period P - and a reliable 
code must be able to stably evolve stable rotating stars 
for several of both timescales. 

Results for star C with Tff — 1.83 and P = 26.38 in 
our units, are plotted in Figures 0and|l^. In Fig. |l^ we 
compare the evolution for maximal slicing (|4^), harmonic 
slicing (^H) and the K-driver (^l|), all with the Gamma- 
driver shift condition (p9[). We find that there is little 
sensitivity to the lapse choice except for small oscillations 
in J which arc only present with harmonic slicing. 

In Fig. ^ where we compare the frozen shift condition 
(H), the AMD shift @ and the Gamma driver @, aU 
evolved with the K-driver (|4^) for the lapse. This com- 
parison demonstrates the great importance of choosing 
an appropriate shift condition for controlling F*. AMD 
is dramatically better than frozen shift in this regard, 
and the Gamma driver is dramatically better than AMD. 
The behavior with AMD shift does not change signifi- 
cantly when the criteria for convergence of (|4j) is made 
stricter. Note that the modification (E^) to AMD and 



FIG. 12: The evolution of the central rest density and the 
ADM mass as star D collapses. 
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FIG. 13: The evolution of the lapse and conformal exponent 
at the origin as star D collapses. 



the Gamma-driver is not activated for this application. 

Figures 12 and |l^ show the behavior of the radially 
unstable star D under different coordinate choices. Once 



again, perturbations are induced solely by the finite dif- 
ference error of the grid. We terminate simulations when 
mass conservation is violated by 10% or the code crashes. 
The singularity avoidance property of the K-driver, which 
approximates maximal slicing, is manifest: with the lapse 
collapsing to very small values, the proper time between 
time slices at the star's center becomes very small, which 
effectively "freezes" all quantities there. With harmonic 
slicing, a decreases more slowly, and we are able to reach 
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FIG. 14: The central density and angular mornonturn of star 
C evolved on a 64 x 32^ grid with K-driver lapse and Gamma- 
driver shift in the inertial and in the rotating frame. We see 
that J is conserved much better in the rotating frame. 

higher central densities, corresponding to later proper 
times, before the code crashes. Given their qualitatively 
different behavior, it is difficult to compare meaningfully 
the different lapse choices for this scenario. If one wants 
to see the central region reach the farthest stage of col- 
lapse before violation of mass conservation becomes un- 
acceptable, harmonic lapse and Gamma-driver shift seem 
to be the optimal combination. One possible reason for 
this is the behavior of 0c, the conformal exponent at the 
stellar center. For the gauge choices which are best suited 
to probing the central region, 0^ decreases significantly 
from its initial value. Inverting Shibata's reasoning for 
modifying the AMD gauge, we infer that this corresponds 
to choosing a gauge with infalling coordinates. This ef- 
fectively increases the grid coverage of the collapsing star, 
resulting in a more accurate evolution. 

We are only able to locate an apparent horizon in 
the harmonic lapse/Gamma driver simulation, and only 
in the last few timesteps, at which point ac — 0.05, 
Mah/M ~ 0.58, and the error in ADM mass is about 
2%. It seems that 64 x 32^ zones are barely sufficient 
resolution for locating horizons for rotating stars reliably. 



2. Rotating Frame 

We compare results for uniformly rotating star C in the 
inertial frame to results in the corotating frame in Figure 
p^ . In the corotating frame, = at t = 0. The light 
cylinder, where points of fixed coordinate label are mov- 
ing on null paths, is at rcyi — 4.2, well outside our outer 
boundaries at x,y, z = 1.5. All coordinate observers are 
therefore timelike everywhere on our grid. We see a dra- 



FIG. 15: Star E evolved for 4 central periods on a 64 x 32^ 
grid. The M and J curves overlap. 



matic improvement in angular momentum conservation 
in the corotating frame. This indicates that the loss of J 
in the inertial frame is caused by error in the advection of 
fluid quantities along v"^. Mass conservation is also bet- 
ter in the rotating frame, but not dramatically. Other 
quantities show qualitatively the same behavior in both 
frames. We have redone the evolution of collapsing star 
D with harmonic lapse and Gamma-driver shift in the 
corotating frame, and our results were almost identical 
to those in the inertial frame. 



C. Differentially Rotating Stars 

We now test the ability of our code to handle differ- 
ential rotation. Differential rotation in neutron stars is 
relevant in several important astrophysical phenomena. 
Simulations in both Newtonian hydrodynamics p8[ and 
full general relativity pl| , ^ indicate that binary neu- 
tron star coalescence may well lead, at least temporar- 
ily, to a differentially rotating remnant, which can sup- 
port significantly more rest mass than uniformly rotating 
stars |23| . Core collapse in a supernova may also result 
in a differentially rotating neutron star. 

We construct axisymmetric equilibrium initial data, 
again following p5[ |, with z chosen as the axis of sym- 
metry. For the rotation profile, we choose 

u'u^ = Rl^A''{n,~n) (67) 

where il is the angular velocity of the fluid, fie is the 
value of on the rotation axis, i?oq is the equatorial 
coordinate radius, and A is a parameter that measures, 
in units of Rcq, the scale over which n changes. In the 
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FIG. 16: The evolution of star F with 0%, 90%, and 99.9% 
of its pressure removed, respectively. When no pressure is 
removed, the star is stationary. When 90% is removed, we 
evolve until the new equilibrium is reached. When 99.9% is 
removed, the star collapses from an initial radius of 26M to 
a radius of ~ 5A/, at which point the simulation becomes 
inaccurate, and we terminate it. 

Newtonian limit this profile reduces to 

For A — > oo one recovers uniform rotation. 

In Fig. |l^ we present results for star E with Pmax — 
0.07, A-^ = 1, i?oq/M = 4.48, T/\W\ = 0.23, and 
J/M^ = 1.02. This star's rest mass of Mq = 0.304 
exceeds the maximum allowed rest-mass of non-rotating 
r = 2 polytropes by 70 percent. We evolve this star on a 
64 X 32^ grid, using 7r-symmetry, with outer boundaries 
at [—2,2] X [0,2]^. The same star was evolved dynami- 
cally by [ p3[ , and we confirm their finding that the star 
is stable over several central rotation periods. 

We found that simulations of differentially rotating 
stars are very sensitive tests of hydrodynamic advection 
schemes. In particular, when we used time-averaging in- 
stead of Crank-Nicholson to treat the advection terms 
(see Appendix ^) , we found that the angular momentum 
is conserved very poorly. The decrease in J also causes 
the central density to rise, and the numerical model to 
drift further and further from the true solution. This 
suggests that for differential rotation, the ability to suc- 
cessfully conserve angular momentum depends strongly 
on the finite difference algorithm used for the hydrody- 
namics. 

In Fig. hq we show results for star F, with pmax = 
0.0174, A^= 3, Req/M = 26.3, T/\W\ ^ 0.0528, and 
J/AP = 0.715. This model is identical to star I in Table 
II of p9|, where this star was evolved in axisymmetry. 



FIG. 17: The evolution of star F with 99.9% of its pressure 
removed. This time, we evolve on a 100^ X 50 grid. The points 
mark times when the resolution was doubled. 



Star F is radially stable, but, as in we make the 
situation dynamic by depleting pressure from the star 
by artificially reducing the polytropic constant k (which 
requires us to re-solve the Hamiltonian and momentum 
constraints). Removing pressure support causes the star 
to implode. For small depletion factors, this collapse is 
halted and the star bounces and finds a new, more com- 
pact, stable equilibrium configuration. When k is re- 
duced to a low enough value, the star will collapse to a 
black hole. Thus, there is a critical polytropic constant 
Kcrit separating these two outcomes. In [ jSOf , this critical 
value was found to be Kcrit ~ 0.04. 

We evolve star F on a grid of 64^ x 32 zones, with 
the outer boundaries located at 2, or equivalently 40.8 
M . In this simulation we use only equatorial symmetry 
so that non-TT-symmetric perturbations can grow. We 
evolve three different cases; one without pressure deple- 
tion with K = 1, a supercritical case with k = 0.1 > Kcrit, 
and a subcritical case with k — 0.001 < Kcrit- Both 
the second and third case present unique challenges. In 
the second case, the collapse is halted by a strong shock 
which must be handled accurately. In the third case, 
we must follow the collapse from a radius of 26.3Af to a 
radius of « M. Our results are consistent with those of 
, even though our 3D simulations have a much poorer 
resolution than the axisymmetric simulations of [ p9| . In 
particular, our resolution is insufficient to follow the final 
stages of the k = 0.001 collapse and prove that a black 
hole is formed. 

In order to overcome this problem, we redo the k — 
0.001 collapse on a 100^ x 50 grid. This grid is still too 
sparse to resolve a black hole with radius of approxi- 
mately IM if the outer boundaries are imposed at 40.8 
M . In order to resolve the black hole, we rezone our grid 
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several times during the implosion, halving the bound- 
aries and halving the grid spacing, so that the total num- 
ber of gridpoints remains constant (compare We 
present results for a simulation that was carried out on 
four different grids with outer boundaries at 2, 1 0.5 and 
0.25. We use the K-driver and the Gamma-driver without 
the modification (^). With the modification, the func- 
tions r* grow very rapidly and cause the code to crash 
well before the radius reaches M. Turning off the modifi- 
cation means that 0c will grow, and the coordinates will 
blow outward. We count on the grid rezoning to counter 
this effect. Also, we switch to frozen shift on the last and 
finest grid of the evolution. The Gamma-driver does not 
perform well on this segment, perhaps because the grid 
boundaries have been moved in to a point where the shift 
does not have its asymptotic form (p8|). 

The results of this simulation are shown in Figs, 
and m. M and J remain within 10% of their initial val- 
ues throughout (we terminate the calculation when this 
ceases to be true.) In our coordinates, the equatorial ra- 
dius decreases from 1.24 (25.3M) to 0.04 (0.8M). Since 
(j)c is growing, part of this decrease in radius is a coor- 
dinate effect. The coordinate- independent circumferen- 
tial radius at the equator (computed from g^^) decreases 

1 /2 

from 27M to 1.7M. At tpc — 0.98, we locate an appar- 
ent horizon with surface area A — 0.0804 corresponding 
to an "irreducible" mass of {A/lGn^Y^'^ = 0.8M. Is this 
a reasonable value? The existence of rotation and of mass 
outside of the black hole mean that we can no longer ex- 
pect the irreducible mass of the hole to be equal to the 
ADM mass of the entire system. The area of the event 
horizon of a Kerr black hole with this system's total M 
and J would he A = 0.109. By breaking up the rest 
mass integral into pieces inside and outside the horizon, 
we find that 82% of the baryonic mass is inside the ap- 
parent horizon. If we assume that the values of M and 
J/M for the black hole are 82% of those of the total sys- 
tem, we arrive at the very crude estimate A — 0.0732, 
which is within 10% of the value determined from the ap- 
parent horizon. We terminate our simulation 2.5M after 
the horizon is located, during which time its area does 
not change appreciably. 

Our agreement with indicates that nonaxisymmet- 
ric perturbations are not important in the collapse of this 
star. We confirm this in Fig. |l|. As one can see, the den- 
sity profiles remain axisymmetric throughout. 

IX. BINARY SYSTEMS 

Binary neutron stars are among the most promising 
sources of gravitational radiation for the new generation 
of gravitational wave interferometers. This makes the 
numerical simulation of such systems one of the most im- 
portant goals of a fully relativistic hydrodynamics code 
and provides one of the most demanding tests for any 
such code. A binary system allows us to uncover poten- 
tial problems that may not be evident in axisymmetric 



scenarios. Previous simulations have focused on the co- 
alescence and merger of binary neutron stars [^l], 
In this Section we demonstrate that our code can stably 
evolve binaries in stable, quasi-circular orbits for over two 
periods (compare [p7|). 

As initial data for these simulations we adopt the data 
of |]6l| , |62| , describing two equal mass polytropes in co- 
rotating, quasi-circular orbit. These data have been con- 
structed using the conformal "thin-sandwich" decompo- 
sition of the constraint equations ||ll|, |2[ |l^, |6^, Q to- 
gether with maximal slicing and spatial conformal flat- 
ness. 

In this Section we focus on one particular case and 
postpone a more complete presentation for a forthcoming 
paper Q. We model the neutron stars as F = 2 poly- 
tropes with an individual rest mass of Mq"'^ = 0.1 in our 
nondimensional units (recall that the polytropic index k 
is set to unity). At infinite separation, this corresponds 
to an individual gravitational mass of M™'^ = 0.096. The 
compaction of (M™'^/i?)oo — 0.088 implies that the grav- 
itational fields are moderately relativistic (the maximum 
compaction for F = 2 polytropes is (Af'"'^/_R)oo = 0.21). 
We adopt initial data for a binary separation of Za — 0.3, 
where Za is the ratio between the coordinate separation 
from the center of mass to the nearest point on the star's 
surface to the farthest point (see |6^, |6^), meaning that 
the separation between the stellar surfaces is about 86% 
of a stellar diameter. This separation is well outside the 
innermost stable circular orbit (ISCO) as determined by 
the analysis of initial data sets (see |^). At this sepa- 
ration, the total binary ADM mass is M = 0.19 and the 
total angular momentum is J/M"^ = 1.36. 

We evolve these initial data on three different grids. 
Two "small grid" simulations are evolved on 120 x 60^ 
gridpoints, with a resolution of Aa; = Ay — Az — 0.55M 
(the binary is symmetric across the equatorial plane, and, 
for equal mass stars, 7r-rotation symmetric around the 
center of mass). The individual stars are resolved by 
~ 16 gridpoints across the stellar diameter (compare |22f| 
where much larger grids are used). One of these small 
grid evolutions is performed in the inertial frame, the 
other in a rotating frame. On these small (uniform) 
grids, the outer boundaries are imposed very close to the 
stars (at a separation of two stellar diameters), which 
we found to introduce numerical noise. We therefore re- 
peated these simulations on a "large grid" , performed in 
a rotating frame, where we doubled the number of grid- 
points and the separation to the outer boundary, while 
keeping the grid resolution constant. This corresponds 
to a numerical grid covering a cubic coordinate volume 
of side [-66,66] in the units of Fig. |l^. The size of this 
grid is such that the corner points almost "touch" the 
surface of the light cylinder, the cylinder with coordinate 
radius Rl = where is the rotating frame angular 
frequency. The possibility of evolution in rotating frames 
with grids that extend beyond the light cylinder will be 
studied in a future article. 

We used Courant factors of 0.30 and 0.46 for the small 
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FIG. 18: Snapshots of the rest density contour hnes for po and the velocity field (t;^, u^) in the equatorial plane for a simulation 
of the collapse of star F with k. = 0.001. The contour hnes are drawn for po = 10^<°-^ 3+o.i)pMax^ where po'"'' denotes the 
instantaneous maximum value of po for j — 0, 1, .., 7. Vectors indicate the local velocity field, v^. The thick circle on the last 
frame marks the apparent horizon. 



and large grid runs respectively, resulting in about 3,000 
timesteps per orbit for the latter case. We adopt the 
K-driver ([4l| ) for the lapse (using e = 0.625, c = 0.1, 
and 5 substeps per step) and the gamma-driver ( ^9| ) for 
the shift (using rj = 0.2, k — 0.005 and 10 substeps). 
We used Sommerfeld type bou ndary conditions for the 
gravitational fields (see Section [ID) an d Cop y type for 
the hydrodynamical fields (see Section HIE). For the 
artificia l visc osity we used CqvIs =0.1 and ClvIs = (see 
Section pT^ ). We also used ch = 0.050Ar in (|l|). For 
the small, rotating frame grid, this choice led to the code 
crashing after about a period and a half. However, we 
found that restarting the code just before that with ch = 
0.024AT allowed us to continue the evolution. For the 
large grid, no such adjustment was necessary. For both 
small grid simulations we found that the stars quickly 
drift apart. To compensate for this we reduced the orbital 
angular velocity fl by 2% for these two cases. Again, no 
such adjustment was required for the large grid. 

In Fig. ^ we show contour plots of the rest density 
Po at half period intervals for the large grid simulation 
in the rotating frame. The three-velocity of the fluid 
is represented by the arrows. The fact that the different 
panels look almost identical indicates how well the binary 
remains in its circular orbit. 



Imposing the outer boundaries at smaller separations 
we were unable to keep the binary in circular orbit. In 
Fig. 20 we plot the coordinate separation d between the 
two points of maximum density for all three simulations 
as a function of time [|67| . For the small grid in the iner- 
tial frame, the two stars start to drift apart after a very 
short time. When performing the same simulation in a 
rotating frame, the stars remained in binary orbit for 
about two periods, but ultimately merge. This merger 
is triggered by the development of orbital eccentricity. 
When we compare this small grid run with the large grid 
simulation, we see that the eccentricity is greatly reduced 
|]68| . This result seems to validate the quasi-equilibrium 
approach to obtaining reasonable initial data for coro- 
tating neutron star binaries in circular equilibrium and 
underlies the importance of the boundary proximity ef- 
fect in these simulations. 

We show more diagnostics of these runs in Figs. |2l| 
through 25, In Fig. we show the rest mass Mq, gravi- 
tational mass M, and angular momentum J for the three 
different simulations. The use of a rotating frame, which 
minimizes fluid advection through the numerical grid, 
leads to large improvements, especially in the conserva- 
tion of angular momentum. Close outer boundaries lead 
to considerable numerical error. The system loses mass 
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FIG. 19: Snapshots in a rotating coordinate frame of the rest density contour lines and the velocity field in the equatorial 
plane for a simulation of a corotating binary. The contour lines are drawn for po — 10^'"'^ "'^"'^Vof"'^: where Pof""^ denotes the 
maximum value of the rest-density po at t = (here it is 0.0573), for j = 0, 1, .., 7. Vectors indicate the local velocity field and 
the scale is as shown in the top left-hand frame. The stars are orbiting clockwise in the inertial reference frame with an initial 
coordinate velocity of 0.102c. P denotes the initial orbital period. 
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FIG. 20: Evolution of the coordinate separation between 
the maximum rest mass density points of the two stars in 
the binary system shown in Fig. FL3. The time is given as a 
fraction of the initial orbital period and the separation d as a 
fraction of the initial value di. 
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FIG. 21: Rest mass, gravitational mass, and arigular momen- 
tum for the three simulations depicted in Fig. GG . 



and angular momentum through the emission of gravita- 
tional radiation, but at a rate that should lead to smaller 
deviations than we find in our simulations. The maxi- 
mum variation of the rest mass, gravitational mass, and 
angular momentum for the large grid run over the first 
two orbits was 0.3%, 0.3%, and 2.2% respectively. 

In Fig. ^ we show the maximum rest density po 
and the minimum value of the lapse a as functions of 
time. The small oscillations correspond to the funda- 
mental mode of the individual neutron stars, which are 
induced by the truncation error of the finite grid resolu- 
tion. The period of the oscillations P ^ 16.6 agrees well 



FIG. 22: Maximum rest mass density po (top) and minimum 
lapse function a (bottom) for the large grid, rotating frame 
simulation. 
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FIG. 23: L2 norms of the Hamiltonian constraint Ti (top) and 
momentum constraints (bottom) for the large grid, rotat- 
ing frame simulation. All the curves have been normalized as 
explained in Sec. pO. 



with the theoretical value oit = 16.0 (see Fig. 32 in |_7|). 
These oscillations are not damped, since for these runs 
we switched off the linear viscosity terms (ClvIs = 0). 

We monitor the Hamiltonian (p^), momentum (p^), 
and Gamma (|l^ ) constraints in Figs. p3| - ^ . We show 
the L2 norms of the corresponding constraint violations. 
For the Hamiltonian and momentum constraint, these 
violations are normalized with respect to Nhc £^nd N]\jc 
evaluated at / = (see Eqs. (59) and (60)). In Figure 
we show the small-grid rotating-frame result for the 
Hamiltonian constraint violation (solid line), as well as 
the result from a similar small-grid evolution in which 
we set Ch = (dashed line) |6^]. The difference between 
these two lines illustrates the effect of the addition of this 
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FIG. 24: L2 norms of the Hamiltonian constraint H, for small 
grid, rotating frame simulations, using ch = 0.05AT (solid 
curve) and ch = 0.00 (dashed curve). This plot shows the 
effect of the ch term in Eq. ( p^ ) on the conservation of the 
Hamiltonian constraint. 



"to 



0.001 



0.0008 



0.0006 



0.0004 



0.0002 



X comp. 




y comp. 




z comp. 






/ / 

/ ^ 











0.5 



1 

t/P 



1.5 



FIG. 25: L2 norm of the gamma constraint for the large 
grid, rotating frame simulation. 



particular term in Eq. (|l9| ) on the conservation of the 
Hamiltonian constraint. The constraint deviations tj% 
which are particularly sensitive to the choice of spatial 
gauge, remain well behaved during the first two orbits in 
these evolutions with the Gamma-driver. 

As this example demonstrates, our code is able to sta- 
bly evolve binaries in stable, quasi-circular orbits for over 
two orbital periods. In a forthcoming paper |30 we will 
use this code to systematically study binary sequences, 
both to dynamically locate the ISCO and to test how ac- 
curately currently available quasi-equilibrium initial data 
represent binaries in quasi-circular orbit. 



X. SUMMARY 

We have tested our 3-1-1 relativistic hydrodynamics 
code on a variety of problems. We find that our cur- 
rent algorithm, supplemented by driver gauge conditions, 
is rather robust. The grid resources required for stable 
evolution and reasonable accuracy are modest. We accu- 
rately evolve shock tubes, spherical dust collapse, and rel- 
ativistic spherical polytropes. We also evolved uniformly 
and differentially rotating equilibrium polytropes, and 
maintained stable configurations stationary for several 
rotational periods. Two applications carried out with 
our code are particularly significant. First, we examined 
the collapse from large radius of a star with significant 
spin to a Kerr black hole. Second, we evolved stable bi- 
nary neutron stars for several orbits, maintaining quasi- 
circular equilibrium. The first application indicates that 
we can study the effects of angular momentum on gravi- 
tational collapse and on the resulting waveform, an effort 
already initiated in [^9[ . The second application indi- 
cates that we can identify and evolve dynamically stable 
quasi-circular neutron star binaries. This ability can be 
used to locate the ISCO dynamically and to follow the 
transition from an inspiral to a plunge trajectory. In ad- 
dition, dynamic simulation allows us to improve binary 
initial data, for example by allowing initial "junk" grav- 
itational radiation to propagate away. We also hope to 
compute detailed gravitational waveforms form these bi- 
naries, refining the wavetrains reported in [^9[ Q . 

We note that several challenges remain to be addressed 
before there exists a code capable of modeling all the 
gravitational wave sources of current astrophysical inter- 
est. One problem is the need to maintain adequate grid 
coverage of the collapsing star or inspiralling binary while 
still keeping the outer boundaries sufficiently distant, i.e. 
the problem of dynamic range. Adaptive mesh tech- 
niques far more sophisticated than the crude rezoning 
used here may be necessary. A related problem concerns 
gravitational wave extraction, as it currently is not pos- 
sible to place outer boundaries in the wave zone. Finally, 
the formation of black hole singularities in hydrodynamic 
collapse scenarios remains an additional challenge to de- 
termining the late-term behavior of such systems. Special 
singularity-handling techniques, such as excision, need to 
be developed further. 
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APPENDIX A: TREATMENT OF ADVECTIVE 
TERMS IN THE HYDRODYNAMIC EQUATIONS 



requires solving equations of the 
d{qv) 



Solving Eqs. ( p6| 
form 



dt 



dx 



S 



(Al) 



Let us finite difference this equation. Let be the value 
of q at gridpoint i on time level n. Let Ax be the coor- 
dinate distance between neighboring gridpoints, and AT 
be the timestep. The Courant factor is C = AT/ Ax. 
Then we difference (|AlD, say for the predictor step, as 

qr'=€+m{vU/2<lU/2-v'l+l/2€+l/2)+Sn (A2) 

where 

„n n I n n 



"i+\l2 
1i+l/2 



q," + AxVfq/2 if vf > 
gr+i - AxV?+ig/2 if < < 



(A3) 
(A4) 



In (A4), 



where A 



^tl/2<l+^:+l/2l 



otherwise 



i+l/2 



q 



i+l 



qf)/Ax. 



(A5) 

In many van Leer 



schemes, the term q + AxVq/2 in Eq. (A4) is replaced 
by the time-averaging expression q + (Ax — wAT)Vq/2. 
Since we use a predictor-corrector method, we don't need 
to time average. 



APPENDIX B: NEWTONIAN LIMIT OF THE 
EULER EQUATION IN A ROTATING FRAME 

Here we recover the Newtonian Euler equation in a 
rotating frame, as given for instance in ||70|| , taking the 
weak-field limit of the general relativistic Euler equation. 
All variables and differential operators are given in the 
rotating frame, with the exception of "barred" quantities 
that reside in the inertial frame. For simplicity, we will 
work with the general relativistic Euler equation written 
as a function of the variable Uk, defined as 



Uk = huk^%,u%e^^{v' +13') . 

Accordingly, we have 

dt{pi,Uk) + di{pi,Ukv') = 

—ae^'^dkP — PirvPahdka — Pi,Ujdk!3^ 



(Bl) 



-2 T. Ofe0 - 



2wO/i 



dkf , (B2) 



where and h are defined in Sec. (Ill A). Using the con- 
tinuity equation (^) we can re-write the Ihs of equation 
(H) as 



dtipi,Uk) + di{p^Ukv') = Pi,idtUk + v'diUk) ■ (B3) 

To obtain the Newtonian limit, we expand the different 
terms in equation (B2) to first order in the Newtonian 
potential 4>n and the square of the fluid velocity v. For 
instance: 



(3' =P' + inxr) 



Iki - 


Ski 




u° - 


14 


2 


h - 


^ 1 




- 


1 - 




a - 


14 




r)' - 


{n 


X f)^ 



(B4) 



where P' is the shift vector in a inertial frame, f2 is the 
angular velocity of the rotating frame with respect to the 
inertial fra me, and r = (x, y, z) is the position vector. 
The limits (B4) combined with equation (Bl) give 



Uk 



(B5) 



We proceed now to take the Newtonian limit of the rhs 
of Eq. (B2). To do so, we note that 



dtuk -> dtv'' 4 dt{(i X ff = dtv'' 
since dt(l — dtf — 0, and 

v'd.iik w'a, (w*^ 4 (f1 X r)*^ 



These conditions together with equation give 
the Newtonian limit of the Ihs of equation (B2) 



(B6) 



5t(p*Ufc) + 9i(p*Ufcv') 
p* (dtv 



(B7) 
give 

(B8) 



To work on the rhs of equation (B2) to Newtonian 
order, we derive the following limits using Eqs. (B4), 
again keeping only the first order terms in (pN and v'': 



ae'^^dkP 
Pi,vP hadkd 



2p^h 



((u0a)2 - 1) 



p*e 
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-Uiiijdk^''- 



2u°h 

The final term is composed using 



dkP 
P*dk4>N 



. 

I): 



(B9) 



p^iijdkP' ^ p. iv^ + (O X r)^ ]dk{nx r)^ . (BIO) 
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We rewrite the first term above as 



Pi,v^dk{^ X fy 



and the second term as 

p^{n X r)^dkin X fy = -p^[nx{nxf) 



(Bll) 



(B12) 



Combining Q, (|B|), ( p31l| ), and ( |B12D , yields 

(9tf + + (f1 X w)*^) = 

-dkP ~ p.dkcj^N ~ P*{^ X v)'' (B13) 

— X (fl X f)j 

Rearranging terms and replacing by its limit the mass 
density p yields the Ne wto nian limit of the general rela- 
tivistic Euler equation (B2) : 



dtv'' +v'div'' = 



--dkP - dk(j)N - 2(f1 X vf f)^^ ,(B14) 

where the last two terms of the rhs correspond to the 
familiar Coriolis and centrifugal force terms. 



APPENDIX C: ADM MASS AND ANGULAR 
MOMENTUM IN ROTATING FRAMES 

In this Appendix, the "barred" fields represent vari- 
ables in the inertial frame, while the non-barred ones 
are quantities in rotating frames. In Sec. we defined 



the total mass and angular momentum of an asymptot- 
ically flat spacetime by two surface integrals (Eqs. (|6l| ) 
and ( |6^ respectively) which characterize the asymptotic 
behavior of the metric on a time slice. These surface in- 
tegrals were transformed into the volume integrals (163) 
and (Q) according to the calculation described in [pSf . 
These volume integrals are numerically evaluated in our 
code on the computational grid. When working in rotat- 
ing frames, one might worry that these integrals do not 
apply, since the 4-metric is not asymptotically flat due to 
the Vl X r term in the shift. It turns out that this is not a 
problem, since the surface integral formulae for M and J 
can be obtained assuming only that the 3-metric and ex- 
trinsic curvature are asymptotically flat Therefore, 
we can evaluate the volume integrals (|6^) and in the 
rotating frame and be sure that the M and J that we find 
at a given time will be the same as what we would have 
found by transforming into an inertial frame and then 
computing the integrals. We can see this explicitly by 
transforming the integrands from an inertial to a rotat- 
ing frame. For example, the mass (|6^ ) written in terms 
of the "barred" inertial frame quantities is 



M ^ 



IM'^-^^WM'-^.^'^ (CD 



1 ~ 

16^ 



-L i jik 



167r 



For simplicity, we take the inertial coordinate system to 
be the one which is instantaneously aligned with our ro- 
tating frame at the time that we are computing M and 
J. Then the transformation is given by Eqs. ([Hi]), ([s^), 
and (^. (From Eq. (|^), we see that p is an invariant.) 
Applying these rules, we see that every term in the in- 
tegrand is identical in the inertial and rotating frames. 
The same is true of the integrands for J and Mq. 
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